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Abstract. We present the results of numerical integrations yielding the structure 
of and meridional flow in axisymmetric thin viscous accretion disk models. The 
solutions are obtained by simplifying and approximating first the equations, using 
systematic asymptotic expansions in the small parameter e, measuring the rela- 
tive disk thickness. The vertical structure is solved including radiative transfer in 
the diffusion approximation. Carrying out the expansion to second order in e we 
obtain, for low enough values of the viscosity parameter a, solutions containing 



backflows. These solutions are similar to the results first found by Urpin (1984), 



who used approximations that are only valid for large radii and the asymptotic 



analytical solutions of Kluzniak & Kita ( 1997 ) , valid only for polytropic disks. Our 



results may be important for several outstanding issues in accretion disk theory. 
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1. Introduction 

The Shakura-Sunyaev (hereafter SS) model of geometrically thin steady accretion disks 



has extensively been used, since its publication (Shakura & Sunyaev 1972), to model 
disks in a variety of astronomical objects. In particular, accretion disks, widely accepted 
to be present in mass transferring close binary systems (like CV - cataclysmic variables) 
and in young, still accreting, stars (like classical T-Tauri), were modeled using the SS 
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paradigm and the calculated emerging spectrum was compared with observations. The SS 
model, which is based on the now famous a viscosity prescription, was found to be quite 
successful in reproducing at least the gross features of the observations and because of its 
relative simplicity and successful applicability it has achieved the status of a text-book 



paradigm (see e.g. Frank, King & Raine 2002). 

Two key ingredients facilitate the SS approach - the a viscosity prescription and 
the assumption that the vertical extension of the disk is very small as compared to any 
significant radial scale. The second of these is viable if the cooling of the disk material (by 
radiative processes) is very efficient, so that the sound speed Cg is very much smaller than 
the corresponding swirling flow velocity ri7, which itself is close to the Keplerian velocity. 
Choosing a representative typical value for the radius f and denoting the sound speed 
at that radius by c^, we obtain a natural small parameter for the problem e = Cs/(r2f), 
where SI is the Keplerian angular velocity at the representative point. 

The SS model can naturally be understood in this context as the lowest order (in e) 
approximation (see Regev 1983), but it requires also some consideration of the vertical 
structure, because the radiation (at least in the optically thick case) escapes essentially 
in the z direction. The common practice to treat accretion disks using the SS paradigm 
has thus always been accompanied by vertical "averaging" or "integration", that is, an 
effective elimination of the z dependence from the equations. This gave a set of ordinary 
differential equations (ODE) in r whose solutions provide the disk structure. 

Studies wishing to model more accurately phenomena beyond the SS model, like time- 
dependent behavior, radiatively driven mass loss, various instabilities etc. have stretched 
the SS model to its limits, but the trend was to stay within the vertical averaging practice, 
because of the simplifications it introduces. For a fairly recent review of such applications 



see e.g. Lin & Papaloizou ( 1996 ). 

The question of the origin of the viscosity and the correct implementation thereof has 
also been the subject of extensive research efforts. The present day concensus identifies the 



magneto-rotational instability (Balbus & Hawley 1991) as the source of MHD turbulence, 
which ultimately produces angular momentum transport (see Balbus & Hawley 1998 for 
a review). Although it seems that the classical a prescription is consistent with this 



process for vertically averaged thin accretion disks (Balbus & Papaloizou 199E ) , it is still 
unclear what the adequate model is to use in multi-dimensional time-dependent studies. 



as it was recently pointed out by Hawley, Balbus & Stone (2001). This work actually 
addresses nonradiative accretion flows (NRAF), that is, flows in which the radiative 
cooling time is very long compared to the time scales of interest. Previous studies in this 
context introduced the so-called ADAF (Narayan & Yi 1994, Abramowicz et al. 1995) 
and CDAF (see Abramowicz & Igumenshchev 2001 and references therein) giving rise to 
some controversy among the different researchers. 
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Our work here deals with a thin, efficiently radiating, disk and therefore seems to 
have no direct relation with the NRAF and associated problems. Still the robustness of 
our results (which hold in the polytropic and therefore adiabatic case as well) hints that 
backflows may have some consequences also in NRAF. 

In the study of thin disks it is natural to employ systematic asymptotic expansions 



in the small parameter e. Regev (1983), hereafter R83, introduced it for the purpose 
of calculating the structure of an accretion disk-star boundary layers and carried it out 
only in the lowest order. In that work, matched asymptotic expansions were used, and 
it developed into a full semi-analytic study of boundary layers in disks around T-Tauri 
stars. However, not only were these expansions done only up to the first order, but also 
the vertical structure of the disk (and the boundary layer) was averaged out (using the 
usual vertical integration procedures mentioned above). 



Urpin (1984) was intetested in finding an approximation for a steady accretion disk 
flow, beyond the SS solution, which provides only the average inflow velocity. He was able 
to get the two-dimensional approximate structure of the flow in the r — z plane. In doing 
this he demonstrated that consideration of the vertical gradients of the viscous stress 
tensor leads to a class of solutions in which the flow direction in the midplane of the disk 
may be opposite to that in the the layers around it, up to the disk surface. In Urpin's 
work the viscous terms were not treated consistently (neglected in one of the equations, 
but retained in the other) and in addition a net zero angular momentum transport was 
assumed. Therefore he obtained backflows for all values of the viscosity parameter a and, 
in any case, his solutions are actually valid only for large radii. 

Already the first attempts at a full two-dimensional hydrodynamical numerical simu- 
lations of radiative viscous accretion disks demonstrated that indeed the meridional flow 
cannot often be well described by the vertically averaged value of the velocity. In such 
calculations the result obviously depends on the outer boundary conditions. Kley & Lin 



(1992) found backflows near the midplane, similar to Urpin's solutions, but have shown 
that if an inflow condition at the outer boundary is assumed, a circulation pattern is ob- 



tained. More detailed studies by Rozyczka et al. (1994) also exhibited flows with regions 
of midplane backflows and with circulation patterns. 



Kluzniak & Kita ( 1997 ), hereafter KK, in a beautiful work, carried out expansions 
in e up to second order and considered the full viscous flow dynamics, but neglected 
thermal effects (which Urpin retained). Using a polytropic relation and an inner zero 
torque boundary condition, they obtained global analytical solutions with backflows. In 
these solutions the location of the flow stagnation point on the equator Tgtag was found 
as a function of the viscosity parameter a. 

In the present paper we would like to present the results of asymptotic semi-analytical 
calculations of steady thin accretion disks, similar to the ones obtained by KK, but with 
thermal effects explicitly included via an energy equation, containing viscous dissipation 
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and radiative losses (treated in the diffusion approximation) . The inclusion of the energy 
equation will force us to use some (quite straightforward) numerical ODE integrations 
and thus the solutions will not have the closed analytical form of KK. 



2. Assumptions and equations 

We use the standard assumptions of steady state, axial and equatorial symmetry, sim- 
plifying significantly the equations of viscous fluid dynamics when they are expressed in 
cylindrical coordinates. 

Writing the equations in a non-dimensional form, that is, scaling all the variables 
by their typical values (denoted by tilde), brings out the non-dimensional constant e. 
This constant was defined above and it reflects the disk thickness. It is thus a natural 
expansion parameter for a perturbative treatment of geometrically thin (e <C 1) disks. 

Denoting the radial, vertical and azimuthal velocity (in cylindrical coordinates r, z, (p) 
by u,v and rO respectively, the statements of conservation of mass, angular momentum, 
horizontal and vertical linear momentum and energy, in steady state and with the as- 
sumption of axial symmetry, read: 

e - dr{rpu) + dz{pv) = 0, (1) 
r 

where p is the mass density; 

f II 1 

drir^n) + vd,n= -^drivr^ drQ) + -d,{rjd,n), (2) 

where rj = pv is the viscosity coefiicient. We note in passim that this is, of course, the 
anomalous viscosity and it is assumed that irrespective of its source it can be parametrized 
using the usual (SS) a prescription, which allows us to express the viscosity coefiicient 
using the well known Ansatz for the r — (p component of the viscous stress tensor 

\Tr,4>\ = ap, (3) 
where p is the pressure; 

udrU + evdzU- O^r = f^^^ -e^^drP + e /^j^, (4) 

eudrV + v d^v = f^,^^ --^d^p + /4; (5) 



and 

r 1 1 

= 9vis-9"' (6) 



^-{eudrT + vd,T) + T 



e- dr{ru) + dzV 



7- 

where T is the temperature and 7 the adiabatic constant and where we have assumed 
that the disk material obeys the perfect gas equation of state. 

Some additional symbols appear in the above equations and they are: 
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1. The gravitational force (per unit mass) components, given by: 

•'grav j,2\ ^^2) ' •'grav ^ Jgrav 

2. The viscous force (per unit mass) components /^j^; /^j^ and the viscous energy dis- 
sipation rate (per unit mass) q^j^. Since the expressions defining them are somewhat 
lengthy we shall refrain from giving them explicitly here and refer the reader to e.g. 
Tassoul (1978). 

3. The energy loss rate (per unit mass) g~, whose form depends on the physical condi- 
tions. Assuming only radiative losses we get 



p 



r 



(7) 



where A is a non-dimensional constant and F^ , F^ are the appropriate radiative 
flux components. If the disk is optically thick (as we shall assume in this work), the 
diffusion approximation for radiative transfer is viable and the radiative fluxes can be 
expressed in terms of the structure functions, temperature gradient and the opacity 
in the usual way (see R83). 

3. Lowest order solutions - the SS disk 

Expanding all the dependent variables in power series in e and denoting the succesive 
terms by the subscripts 0, 1, 2 . . . it is quite easy to see (R83, KK) that in the lowest 
order we get 

- dr{rpoUi) + dz{poV2) = 0, (8) 
r 

because uq — vq ^ vi = 0, 

rpouidrir'^no) = dr{r]r^dr^o), (9) 
r!o = r3/2, (10) 
dzPo = -Pozr^, (11) 
d,Fo = A-^f]r^{drnof, (12) 
where Fq = Fq . 

The SS solution can be easily recovered from these equations (see R83), but our aim 
here will be to obtain from them also the vertical structure, needed for the next order 
approximation. Like in the SS approach we vertically integrate equations (||) and (^) and 
get 



rpQUidz = const = —to, (13) 
where the constant m is the nondimensional mass transfer rate, and 
mdrir^no) ^ -drinr^drno), (14) 
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where ^ = J^^rjdz. Using now ( p^ ) we can integrate equation (14) on r and use the 
condition that at some r — the zero torque condition is satisfied, that is, is at its 
extremum there (see KK). The result is 

= Im (l - ./^) . (15) 



r J 

The fact that a Keplerian angular momentum profile does not have an extremum and 
that it has to be joined somehow to a sub-Keplerian star gives rise to the boundary layer 
problem, which we shall avoid here by looking only at radii sufficiently larger than r+. 

To get the equations determining the vertical structure at each r we use equations 
( |ll| ) and ( [l^ ) and substitute in the latter the Keplerian fio (from equation [l^) . This gives 

dzV = -pzr~^, (16) 

d.F^^rir-\ (17) 

where we have dropped the subscripts from the functions, for convenience. Treating 
the radiation in the diffusion approximation we get the additional equation 

d,T = -KpT-^F, (18) 

where k is the (properly non-dimensionalized) opacity, assumed known as a function of 
p and T. 

This set of three equations, supplemented by the equation of state p — pT, is sufficient 
in principle for the determination of the vertical structure at each radial point r, but 
appropriate boundary conditions are obviously also needed. This complicates matters, 
because whereas the lower boundary can be simply set at z = 0, the upper one cannot be 
determined a priori. We chose to employ the Eddington approximation for this boundary 
condition. This approximation requires radiative equilibrium, that is no energy sources 
in the outer layers, and to achieve it we tacitly assume that the viscosity coefficient 
approaches zero for a small enough optical depth. Actually, our form of the viscosity 
parameter is obtained from the relation (^), which for the Keplerian (lowest order) 
Q{r) gives 
2 

r] = -ar^/^p. (19) 
3 

Typically, p decreases with z and therefore so does 77. Thus using the Eddington approx- 
imation should be reasonable, at least in the optically thick long as we ignore 
possible instabilities which may result in winds and coronae (see Shaviv, Wickramasinghe 
& Wehrse 1999). 

Explicitly, if we define the outer (vertical) limit of the disk and hence of our integration 
as Zout (see below), the integration of equation ( p^ from z = — Zout to z = Zout gives 

F(zout) = ^f^r-' = ±.rhr~y\r'/^ - 4% (20) 
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Fig. 1. Disk structure for a typical case. Contours of the density (in units of p) are plotted 
in the r — z plane. Note that due to the different scalings the figure actually presents a 
vertically stretched disk by a factor of e^^ ^ 50 

where the second equality follows from (^). This is the upper boundary condition on 
the flux, but the value of Zout has still to be specified. The lower boundary condition on 
the flux is much simpler as it is just F[z = 0) = 0. 

The Eddington approximation links the radiative flux and the temperature at the 
base of the photosphere (rdim = 2/3) and this point is deflned to be at 2; = Zout, since 
the flux above it is constant. In its non-dimensional form, using our units, the Eddington 
condition reads 

r^(2out) = ^F{Zont), (21) 

where f = kph is our optical depth unit. 

The three equations ( |l6[]T^ ) are thus solved as a two-point boundary problem subject 
to the four conditions = at z = 0; (|o|), ( |2l] ) and p = at z = Zout, where the value 
of Zout itself is also determined in an iterative way. 

In figure ^ we show the result for p{r, z) in a typical case [a = 0.2 and electron 
scattering, i.e. constant, opacity case). In this, and all subsequent figures the radial 
coordinate is scaled by f = r_|. and the vertical one by h = er^. 

Figure ^ is the vertical temperature profile T{z) for a representative radius value and 
two values of a: 0.4 (the solid line) and 0.7 (the dotted line). As we shall see in the next 
section, these two cases differ radically in their meridional flow pattern, but we should 
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Fig. 2. Vertical temperature profiles for a typical radius in the disk. The non-dimensional 
temperature is shown as a function of z (in units of er+) for a — 0.4 (solid line) and 
a = 0.7 (dotted line) 

stress that the temperature profiles described here are calculated in the lowest order and 
are, of course, not infiuenced by the fiow. 

The procedure outlined in this section provides the lowest order values of the disk 
structure functions po(^j ^)j Po(?'j ^) a-nd so on, which will be needed to solve for the fiow 
pattern in the next order. 



4. Higher order solutions - the meridional flow pattern 

As explained above, the lowest order contributions to the meridional velocity com- 
ponenets are of order e (i.e. u\) and (i.e. v^) and it turns out that the lowest or- 
der correction to the Keplerian angular velocity is of order (^1-2) - see KK and Urpin 
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( |1984D . The relevant equations, following from the expansions and originating from the 
mass, angular momentum and horizontal momentum conservation equations are 

dr{rpoUi) + rd^{poV2) = 0, (22) 
rpou^drir^no) = dr-ivr'^drno) + r^'d.iT^d.n^), (23) 



3^2 

2rponoi^2 = tttPo - drPo + dz{T]dzUi). (24) 
2r^ 



Despite their appearance these equations can be viewed as ODE in z for any given r, 
since pq, po and flo, as well as 77 are already known functions, from the lowest order solu- 
tion for the vertical structure (as explained in the previous section). In addition, equation 
( ^ ) can be decoupled from the other two, because one of the variables, V2, appears only 
in it. Since equations ( ^3[^ ) are of second order each, four boundary conditions are 
required. Two of them are obvious and are the result of equatorial symmetry: 

d^ui = and S^f^s = 0; at z = 0, (25) 

for any r. The other two conditions result from a regularity requirement, that is, that 
ui and remain bounded for z 00. In KK the density achieved a true zero for a 
finite value of z, owing this property to the polytropic assumption. In our case po decays 
with z and we have found that physically meaningful results are to be expected only if 
we require that both poUi and po^^2 are sufficiently small at Zout- This expectation has 
been strengthened by testing our algorithm on the polytropic case and getting very good 
agreement with the analytical results of KK. 

After the two equations ( p3[]2^ ) are solved, V2 follows from the solution of the first 
order equation ( p2|) with the initial condition t;2 = at z = 0. 

The result of a representative calculation {a = 0.4 and electron scattering opacity) 
is shown in Figure ^. For this case we obtain a meridional flow pattern with a backflow 
around the equatorial plane. The stagnation radius is in this case rgtag ~ 7.5r+. 

For other values of a the trend is similar to the one found by KK, that is, rgt^g 
growing with a and for a > acrit there are no solutions with backflows (rgtag — 00). 
While the disk structure depends somewhat on the opacity law (we have experimented 
with electron scattering opacity and Kramers opacity law), the meridional flow pattern 
and the stagnation radius (see below) are quite insensitive to it. In fact the relevant 
graphs (see below) for the abovementioned two opacity laws look indistinguishable and 
consequently we shall display only the electron scattering opacity case. 

A case in which the stagnation radius is so large that it is not found up to r lOOOr^ 
is shown in fig. ^, where a = 0.7. Note that even for these large values of r the flow is 
directed inward everywhere. 

The determination of the value of acrit can be achieved by calculating rstag for a 
sequence of a values and finding where 'rstag(a) grows beyond any reasonable bound. 
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Fig. 3. Meridional flow pattern with backflow. The flow is represented by the direction 
of streamlines at a grid of points in the r — z plane. Note that the scale of z is "stretched" 
by a factor e^^ ~ 50 

Reasonable, in this context, means a radius beyond which the disk does not exist in a 
realistic astrophysical setting. In any case, when r is very large the present numerical 
calculation cannot give reliable results because of the truncation error. 

In Fig. ^ we show the stagnation radius for a number of a values and clearly see 
the trend that it grows beyond any reasonable bound for a <; 0.7. It is interesting to 
note that our results are completely consistent with the analytical result of KK (solid 
line) found for poly tropes. The meridional flow behavior seems to result thus from the 
dynamics rather than from thermal processes in the disk. 

5. Summary 



In this work we have included thermal effects in the scheme, previously proposed by 
Kluzniak and Kita (1997), to calculate the vertical structure of thin accretion disks (in- 
cluding the flow pattern) by perturbative analysis. The small parameter of the expansions 
is e, reflecting the disk thickness. To the lowest order in e the SS solution is recovered 
and we obtain also the vertical pressure and temperature structure. The accretion flow 
is found from the next order and includes, for small enough values of the viscosity pa- 
rameter, equatorial backflows. 
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Fig. 4. Meridional flow pattern with no backflow. The flow is represented by the direction 
of streamhnes at a grid of points in the r — z plane. Note that the scale of z is "stretched" 
by a factor e^^ ~ 50 

These findings support the claim made by KK that the backflows are of dynamical, 
rather than thermal, origin and probably result from vertical gradients of the viscous 
stress tensor. Our result for the flow is obtained numerically, using two explicit initial 
condition at the disk mid-plane and two regularity conditions. The latter complete the 
constraints needed for the solution of a fourth order differential system. However, because 
these constraints are independent of r, our solutions are obtained for all r values and 
cannot account for given detailed boundary conditions on the flow at some r = rout- 

Still, it is conceivable that if such conditions are imposed in a full numerical hydro 
simulation, our solution will be correct only for an intermediate region, situated between 

and rout, and not too close to these boundaries. Our solutions should smoothly join 
the outer flow pattern, which has to be consistent with the outer boundary conditions. 
For an outer boundary condition of inflow along all z, this should give rise to circulation 
patterns of the type observed in the numerical simulations of radiative disks cited in the 
introduction. 

Although our results are formally valid only for very small e and for radiative (as well 
as polytropic, see KK) disks, it is reasonable (as is often the case in asymptotic analyses 
of this type) that they carry over also to cases of e close to unity. This means that quasi 
spherical flows, like the non-radiative ones (NRAF), may be endowed with circulation 
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Fig. 5. Logarithm of the stagnation radius (in units of r+) as a function of a. The circles 
are our numerical results while the solid line is the KK analytical solution for a polytrope. 



patterns as well. Whether these results are directly relevant to CDAF (if these flows 
are indeed proven to be viable) remains unclear. It can, however, be definitely stated 
that ADAF solutions using vertical averaging may not faithfully represent radial energy 
transport when there are substantial backflows, which naturally carry energy outwards in 
the hottest parts of the disk. In addition, we think that an approach of the kind reported 
here may provide useful tools for the understanding of some complicated features of flow 
patterns and instabilities in accretion disks. 
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